#include <cmath>
using std::ceil;

#include <cstdlib>

#include <fstream>
using std::ifstream;
using std::ios;

#include <iostream>
using std::cout;
using std::endl;

#include <list>
using std::list;

#include <map>
using std::map;

#include <ostream>
using std::ostream;

#include <queue>
using std::queue;

#include <sstream>
using std::stringstream;

#include <string>
using std::basic_string;
using std::char_traits;
using std::string;

#include <utility>
using std::make_pair;
using std::pair;

#include <vector>
using std::vector;


#include <tbb/atomic.h>

#include "tbb/blocked_range.h"
using tbb::blocked_range;

#include "tbb/compat/thread"

#include "tbb/concurrent_hash_map.h"
using tbb::concurrent_hash_map;

#include "tbb/concurrent_vector.h"
using tbb::concurrent_vector;

#include "tbb/mutex.h"
using tbb::mutex;

#include "tbb/parallel_for.h"
using tbb::parallel_for;

#include <tbb/parallel_reduce.h>

#include "tbb/parallel_sort.h"
using tbb::parallel_sort;

#include "tbb/pipeline.h"
using tbb::filter;
using tbb::flow_control;
using tbb::make_filter;
using tbb::parallel_pipeline;

#include "tbb/task_scheduler_init.h"
using tbb::task_scheduler_init;

#include "tbb/tbb_allocator.h"
using tbb::tbb_allocator;


typedef basic_string<char, char_traits<char>, tbb_allocator<char> > String;

typedef string Sequence_name;

typedef unsigned char Letter;
typedef vector<Letter> DNA_sequence;

typedef pair<Sequence_name, DNA_sequence> Named_DNA_sequence;

typedef vector<int> Blocked_sequence;

typedef pair<pair<size_t, size_t>, pair<size_t, size_t> > Match;
typedef vector<Match> Matches;
typedef map<Sequence_name, Matches> Results;




/* Represents symbols of sequence */
struct DNA
{
	static Letter const A;
	static Letter const T;
	static Letter const G;
	static Letter const C;
	static Letter const NO_SYMBOL;
};

Letter const DNA::A = 0;
Letter const DNA::T = 1;
Letter const DNA::G = 2;
Letter const DNA::C = 3;
Letter const DNA::NO_SYMBOL = 4;




/* Represents block of _size characters. */
class Block
{
public:
	static size_t const DEFAULT_SIZE;
	
	/* Creates block with default size */
	Block()
		: _size( 10 )
		,_current_size( 0 )
		, _value( 0 ) 
	{
		set_masks();
	}
	
	/* Creates block with size 'size' */
	Block( size_t size )
		: _size( size )
		, _current_size( 0 )
		, _value( 0 )
	{
		set_masks();
	}
	
	/* Creates block with size 'size' and value 'value' */
	Block( size_t size, size_t value )
		: _size( size )
		, _current_size( size )
		, _value( value )
	{
		set_masks();
	}
	
	
	/* Returns value of the block */
	size_t value() const
	{
		return _value;
	}
	/* Adds the 'symbol' to the end of block and drops first block symbol (if exists) */
	void next_symbol( Letter symbol );
	
	/* Returns _size that is the maximum number of characters that block can contain */
	size_t max_size() const
	{
		return _size;
	}
	
	/* Returns _current_size that is count of symbols that are now in block */
	size_t size() const
	{
		return _current_size;
	}
	
private:
	/* Sets masks for this  */
	void set_masks();
	
private:
	size_t _size;
	size_t _current_size;
	size_t _value;

	size_t _extract_value_mask;
	size_t _change_value_mask;

	static size_t const LETTER_MASK;
	static size_t const SIZE_OF_LETTER;
};

size_t const Block::DEFAULT_SIZE = 10;
size_t const Block::LETTER_MASK = 3;
size_t const Block::SIZE_OF_LETTER = 2;

void Block::next_symbol( Letter symbol )
{
	if( symbol < DNA::NO_SYMBOL )
	{
		if( _current_size < _size )
		{
			++_current_size;
		} 
		_value <<= SIZE_OF_LETTER;
		_value += symbol;
		_value &= _change_value_mask;
	}
}

void Block::set_masks()
{
	_extract_value_mask = 1;
	_change_value_mask = 1;
	
	for( size_t i = 0; i != SIZE_OF_LETTER * ( _size - 1 ); ++i )
	{
		_extract_value_mask <<= 1;
		++_extract_value_mask;
		_change_value_mask <<= 1;
		++_change_value_mask;
	}
	_extract_value_mask <<= 1;
	++_extract_value_mask;
}




/* Reference sequence and test sequences are processed by this class */
class Processor
{
public:
	Processor( size_t threads_num, size_t min_match_len );
	
	/* Sets the given sequence as reference sequence and preprocesses it */
	void set_reference_sequence( DNA_sequence * sequence );
	
	/* Processes given sequense as new input sequence and writes the result of processing to _results array */
	void add_new_input_sequence( Named_DNA_sequence const & sequence );
	
	/* Makes sorting of results */
	void sort_results();
	
	/* Constant getter for _results field */
	Results const & results() const
	{
		return _results;
	}
	
private:
	int _threads_num;
	size_t _min_match_len;
	size_t _block_size;

	Results _results;
	mutex _results_lock;

	DNA_sequence const * _reference_str;
	vector<size_t> _reference_blocked_str;
	vector<vector<size_t> > _reference_hash;
	
	class Reference_processing_functor;
	class Input_processing_functor;
};

Processor::Processor( size_t threads_num, size_t min_match_len )
	: _threads_num( threads_num )
	,_min_match_len( min_match_len )
{
	size_t max_block = ( _threads_num > 12 ) ? 8 : 6;
	_block_size = ( min_match_len < max_block ) ? min_match_len : max_block;
}


class Processor::Reference_processing_functor 
{
public:
	Reference_processing_functor( Processor * processor, vector<mutex> * locks ) 
		: _processor( processor )
		, _locks( locks )
	{}
	
	void operator() ( blocked_range<size_t> const & r ) const
	{
		Block block( _processor->_block_size );
		size_t start = _processor->_block_size - 1 + r.begin();
		for( size_t i = r.begin(); i < start; ++i )
		{
			block.next_symbol( ( *_processor->_reference_str )[i] );
		}
		
		size_t end = r.end() + _processor->_block_size - 1;
		for( size_t i = start; i < end; ++i )
		{
			block.next_symbol( ( *_processor->_reference_str )[i] );
			
			( ( *_locks )[ block.value() % _locks->size() ] ).lock();
			_processor->_reference_hash[ block.value() ].push_back( i - _processor->_block_size + 1 );
			( ( *_locks )[ block.value() % _locks->size() ] ).unlock();
		}
	}
	
private:
	Processor * _processor;
	vector<mutex> * _locks;
};


void Processor::set_reference_sequence( DNA_sequence * sequence )
{
	_reference_str = sequence;
	
	static size_t locks_count = _threads_num * 1024;
	size_t blocks_count = sequence->size() - _block_size + 1;
	size_t size = 1 << ( _block_size * 2 );
	_reference_hash.resize( size );
	vector<mutex> locks( locks_count, mutex() );
	
	size_t interval = static_cast<size_t>( ceil( static_cast<double>( blocks_count ) / static_cast<double>( _threads_num ) ) );
	
	parallel_for( blocked_range<size_t>( 0, blocks_count, interval ), Reference_processing_functor( this, &locks ) );
}

class Processor::Input_processing_functor
{
public:
	Input_processing_functor( Processor const * processor, mutex * m, DNA_sequence const & sequence, Matches * res )
		: _processor( processor )
		, _mutex( m )
		, _sequence( sequence )
		, _res( res )
	{}
	
	void operator() ( blocked_range<size_t> const & r ) const
	{
		Block block( _processor->_block_size );
		size_t start = _processor->_block_size - 1 + r.begin();
		
		for( size_t i = r.begin(); i < start; ++i )
		{
			block.next_symbol( _sequence[ i ] );
		}
		
		size_t end = r.end() + _processor->_block_size - 1;
		for( size_t i = start; i < end; ++i )
		{
			block.next_symbol( _sequence[ i ] );
			size_t idx_seq_start = i - _processor->_block_size + 1;
			
			if( !_processor->_reference_hash[ block.value() ].empty() )
			{
				vector<size_t> const & lst = _processor->_reference_hash[ block.value() ];
				for( vector<size_t>::const_iterator it = lst.begin(); it != lst.end(); ++it )
				{
					if( idx_seq_start != 0 && ( *it ) != 0 && _sequence[ idx_seq_start - 1 ] == ( *_processor->_reference_str )[ *it - 1 ] )
					{
						continue;
					}
					
					size_t idx_ref = *it;
					size_t idx_seq = idx_seq_start;
					
					while( idx_ref < _processor->_reference_str->size()
						&& idx_seq < _sequence.size()
						&& ( *_processor->_reference_str )[ idx_ref ] == _sequence[ idx_seq ] )
					{
						++idx_seq;
						++idx_ref;
					}
					
					--idx_seq;
					--idx_ref;
					
					if( ( idx_seq - idx_seq_start + 1 ) >= _processor->_min_match_len && ( idx_ref - *it + 1) >= _processor->_min_match_len )
					{
						_mutex->lock();
						_res->push_back( make_pair( make_pair( *it + 1, idx_ref + 1 ), make_pair( idx_seq_start + 1, idx_seq + 1 ) ) );
						_mutex->unlock();
					}
				}
			}
		}
	}
	
private:
	Processor const * _processor;
	mutex * _mutex;
	DNA_sequence const & _sequence;
	Matches * _res;
};


void Processor::add_new_input_sequence( Named_DNA_sequence const & sequence ) 
{
	_results_lock.lock();
	Matches & results = _results[ sequence.first ];
	_results_lock.unlock();
	
	size_t blocks_count = sequence.second.size() - _block_size + 1;
	mutex m;
	
	parallel_for(
		blocked_range<size_t>( 0, blocks_count ),
		Input_processing_functor( this, &m, sequence.second, &results )
	);
}




/* Operator << defined for output Processor as a set of results that it contains */
ostream & operator<<( ostream & out, Processor const & processor )
{
	for( Results::const_iterator it = processor.results().begin(); it != processor.results().end(); ++it )
	{
		out << it->first << "\n";
		Matches const & res = it->second;
		for( size_t i = 0; i < res.size(); ++i )
		{
			bool flag = true;
			size_t j = i + 1;
			while( flag && j < res.size() ) 
			{
				if( ( res[ j ].first.second - 1 == res[ i ].first.second
						&& res[ j ].first.first - 1 == res[ i ].first.first
						&& res[ j ].second == res[ i ].second )
					||
					( res[ j ].second.second - 1 == res[ i ].second.second
						&& res[ j ].second.first - 1 == res[ i ].second.first
						&& res[ j ].first == res[ i ].first ) 
					)
				{
					flag = false;
				} 
				else if( res[ j ].first.second - 1 > res[ i ].first.second )
				{
					break;
				}
				++j;
			} 

			if( flag ) 
			{
				out << res[ i ].first.first << " " << res[ i ].first.second << " " 
					<< res[ i ].second.first << " " << res[ i ].second.second << "\n";
			}
		}
	}
	return out;
}

class Results_order
{
public:
	bool operator() ( Match const & left, Match const & right ) const 
	{
		if( left.first.second != right.first.second )
		{
			return left.first.second < right.first.second;
		}
		return left.second.second < right.second.second;
	}
};

void Processor::sort_results()
{
	for( Results::iterator it = _results.begin(); it != _results.end(); ++it )
	{
		parallel_sort( it->second.begin(), it->second.end(), Results_order() );
	}
}




/* Reads one sequence from given ifstream 'stream' to 'buffer' */
Named_DNA_sequence & get_dna_sequence( ifstream & stream, Named_DNA_sequence & buffer )
{
	buffer.first.clear();
	buffer.second.clear();
	Sequence_name & name = buffer.first;
	DNA_sequence & sequence = buffer.second;
	
	if( stream.get() == '>' )
	{
		stream >> name;
	}
	
	while( !stream.eof() )
	{
		int c = stream.get();
		if( c == '>' )
		{
			stream.unget();
			break;
		}
		else
		{
			if( c == 'A' )
			{
				sequence.push_back( DNA::A );
			}
			else if( c == 'T' )
			{
				sequence.push_back( DNA::T );
			}
			else if(c == 'G')
			{
				sequence.push_back( DNA::G );
			}
			else if( c == 'C' )
			{
				sequence.push_back( DNA::C );
			}
		}
	}
	return buffer;
}

/* Processes all input sequences in all given files */
void process_all_test_cases( char ** filenames_start, char ** filenames_end, Processor & processor, int ntoken )
{
	char ** filename = filenames_start;
	ifstream * in = new ifstream( *filename, ios::in );
	
	auto read_f = [&]( flow_control & fc ) -> Named_DNA_sequence * 
	{
		/* If current file's sequences read already close the stream and take next filename */
		if( in->eof() )
		{
			in->close();
			++filename;
		}
		/* If all filenames proceed then stop */
		if( filename == filenames_end )
		{
			fc.stop();
			return NULL;
		}
		/* If took new file for processing then open stream for reading */
		if( !in->is_open() )
		{
			in = new ifstream( *filename, ios::in );
		}
		Named_DNA_sequence * named_dna_sequence = new Named_DNA_sequence;
		get_dna_sequence( *in, *named_dna_sequence );
		return named_dna_sequence;
	};
	
	auto process_f = [&]( Named_DNA_sequence * named_dna_sequence ) -> void 
	{
		processor.add_new_input_sequence( *named_dna_sequence );
		delete named_dna_sequence;
		return;
	};
	
	parallel_pipeline(
		ntoken,
		make_filter< void, Named_DNA_sequence * >( filter::serial_in_order, read_f ) &
		make_filter< Named_DNA_sequence *, void >( filter::parallel, process_f )
	);
	
	delete in;
}


int main( int argc, char * argv[] )
{
	// first command line argument : number of worker threads to use
	int worker_threads_num = atoi( argv[ 1 ] );
	task_scheduler_init init( worker_threads_num );
	
	// second command line argument : window size
	size_t min_match_length = atol( argv[ 2 ] );
	
	Processor processor( worker_threads_num, min_match_length );
	
	// third command line argument : reference sequence file
	ifstream reference_sequence_file( argv[ 3 ], ios::in );
	Named_DNA_sequence reference_dna_sequence;
	get_dna_sequence( reference_sequence_file, reference_dna_sequence );
	processor.set_reference_sequence( &reference_dna_sequence.second );
	
	// following command line arguments : other files containing sequences
	process_all_test_cases( argv + 4, argv + argc, processor, worker_threads_num );
	
	processor.sort_results();
	cout << processor;
	
	return 0;
}
